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A spin-1 Blume-Capel model with dilute and random crystal fields is examined for honeycomb and 
square lattices by introducing an effective-field approximation that takes into account the correla- 
tions between different spins that emerge when expanding the identities. For dilute crystal fields, we 
have given a detailed exploration of the global phase diagrams of the system in ksTd J — D I J plane 
Cy^ . with the second and first order transitions, as well as tricritical points. We have also investigated 

^Nj ' the effect of the random crystal field distribution characterized by two crystal field parameters D j J 

and A/ J on the phase diagrams of the system. The system exhibits clear distinctions in qualitative 
manner with coordination number q for random crystal fields with A/ J, Dj J ^ 0. We have also 
found that, under certain conditions, the system may exhibit a number of interesting and unusual 
phenomena, such as reentrant behavior of first and second order, as well as a double reentrance with 
rsl , three successive phase transitions. 
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I. INTRODUCTION 

Spin-1 Blume-Capel (BC) model^'^ is one of the most extensively studied models in statistical mechanics and con- 
densed matter physics. The model exhibits a variety of multicritical phenomena such as a phase diagram with ordered 
ferromagnetic and disordered paramagnetic phases separated by a transition line that changes from a continuous phase 
transition to a first-order transition at a tricritical point. On the other hand, as an extension of the model, BC model 
with a random crystal field represents the critical behavior of "^He — "^He mixtures in a random media, i.e., aerogel 
where S — and S* = ±1 states represent ^He and *He atoms, respectively^. From the theoretical point of view, 
BC model with a random crystal field (RCF) has been studied by a variety of techniques such as cluster variational 
method (CVM)^, Bethe lattice approximation (BLA)- , effective field theory (EFT)®"^", finite cluster approximation 
(FCA)iiii^, mean field theory (MFT)i^"i^, Monte Carlo (MC) simulationsiS, pair approximation (PA)^, and renor- 
malization group (RG) method^i. Among these studies, EFT and MET have been widely used to investigate the 
thermal and magnetic properties of BC model with a RCF distribution. For example, Kaneyoshi and Mielnicki^ inves- 
tigated the phase diagram of the system for a honeycomb lattice by using EFT with correlations and they found some 
important differences from the results obtained by the standard MET. Similarly, in a recent paper, Yan and Deng — 
considered the same model within the framework of EFT, and they derived the expressions of magnetizations for 



honeycomb and square lattices. On the other hand, in several studies based on MFT ;^'^'^^ the authors paid attention 
to the effects of crystal field dilution on the phase diagrams of the system and they observed that the system may 
exhibit a reentrant behavior, as well as first order phase transitions. 

However, EFT and MFT studies mentioned above have some unsatisfactory results. Namely, the results obtained 
by EFT are limited to second-order phase transitions and tricritical points, and a detailed description of first-order 
transitions has not been reported. Other than this, it is well known that magnetic systems with dilute crystal fields 
exhibit qualitatively similar characteristics when compared to site dilution problem of magnetic atoms. From this 
point of view, for a BC model with diluted crystal fields, MFT predicts that the phase transition temperature of the 
system will remain at a finite value until zero concentration is reached. In this context, we believe that BC model 
with RCF still deserves particular attention for investigating the proper phase diagrams, especially the first-order 
transition lines that include reentrant phase transition regions. Conventional EFT approximations include spin-spin 
correlations resulting from the usage of the Van der Waerden identities, and provide results that are superior to those 
obtained within the traditional MFT. However, these conventional EFT approximations are not sufhcient enough to 
improve the results, due to the usage of a decoupling approximation (DA) that neglects the correlations between 
different spins that emerge when expanding the identities. Therefore, taking these correlations into consideration will 
improve the results of conventional EFT approximations. In order to overcome this point, recently we proposed an 
approximation that takes into account the correlations between different spins in the cluster of a considered lattice^^. 
Namely, an advantage of the approximation method proposed by this study is that no decoupling procedure is used 
for the higher-order correlation functions. 

In this paper, we intent to investigate the effects of RCF distributions on the phase diagrams of spin-1 BC model 
on 2D lattices, namely honeycomb (g = 3) and square (q = 4) lattices. For this purpose, the paper is organized as 
follows: In Sec. |TT] we briefly present the formulations. The results and discussions are presented in Sec. IIIIl and 
finally Sec. IIVI contains our conclusions. 

II. FORMULATION 

In this section, we give the formulation of the present study for a 2D lattice which has N identical spins arranged. 
We define a cluster on the lattice which consists of a central spin labeled 5*0, and q perimeter spins being the nearest 
neighbors of the central spin. The cluster consists of {q + 1) spins being independent from the spin operator S. The 
nearest-neighbor spins are in an effective field produced by the outer spins, which can be determined by the condition 
that the thermal average of the central spin is equal to that of its nearest-neighbor spins. The Hamiltonian describing 
our model is 

i/ = -J^5f5|-^A(5ff, (1) 

where the first term is a summation over the nearest- neighbor spins with S^ — ±1,0 and the term Di on the second 
summation represents a random crystal field, distributed according to a given probability distribution. In this paper, 
we primarily deal with two kinds of probability distributions, namely, a quenched diluted crystal field distribution 
and a double peaked delta distribution which are given by Eqs. ([2|) and ([3]), respectively as follows 

F(A) = p6iD, -D) + {1- p)6iD,), (2) 

P(A) = \{6[D,-{D-A)]+S[D,-{D + A)]}. (3) 

where p denotes the concentration of the spins on the lattice which are influenced by a crystal field D. 

We can construct the mathematical background of our model by using the approximated spin correlation identities^ 
by taking into account random configurational averages 

where /? = l/fc^T, {/.;} is an arbitrary function which is independent of the spin variable Si and the inner (...) and 
the outer (...)r brackets represents the thermal and random configurational averages, respectively. 



In order to apply the differential operator technique^ii^^, we should separate the Hamiltonian ([T]) into two parts as 
H = Hi + H . Here, the effective Hamiltonian Hi includes all the contributions associated with the site i, and the 
other part H does not depend on the site i. 

- H, = ESt + D, {Slf , (6) 

where E — Jj^j ^j i^ ^^^ local field on the site i. If we use the matrix representations of the operators S'f and {S^)"^ 
for the spin-1 system then we can obtain the matrix form of Eq. © 

f E + D \ 

-H,= \ \. (7) 

\ Q -E + D I 

Hereafter, we apply the differential operator technique in Eqs. (j4]) and ([5]) with [fi] ~ 1. From Eq. Q we obtain 
the following spin identity with thermal and configurational averages of a central spin for a lattice with a coordination 
number q as 

((^o)). = //n[l + ^|sinh(JV) + (5|f{cosh(JV)-l}]\\ F(a:)U=o. (8) 

The function F{x) in Eq. ^ is defined by 

F{x)^ I dD,P{Di)f{x,Di), (9) 



where 



/(x,A) = -3 -— ^(^„|5f|^„)exp(/3A„), (10) 

i;„=iexp(^A„)„^^ 

2 sinh(/3x) 



2cosh(/3a;) +e-PDr 

In Eq. (jlOp . A„ denotes the eigenvalues of —Hi matrix in Eq. ([7]), and (p„ represents the eigenvectors corresponding 
to the eigenvalues A^ of —Hi matrix. With the help of Eq. (fTO|) . and by using the distribution functions defined in 
Eqs. ([U and ([3]), the function F{x) in Eq. ^ can be easily calculated by numerical integration. Hereafter, we will 
focus our attention on the construction of the correlation functions, as well as magnetization and quadrupole moment 
identities of a honeycomb lattice with q — i. K brief formulation of the fundamental spin identities for a square lattice 
with q = A can be found in Appendix [^ 

By expanding the right-hand side of Eq. ([S]) for a honeycomb lattice with q = 3, we get the longitudinal magneti- 
zation as 

W. = (('S'o))r = lo + 2,ki{{Si))r+2,{h-lo){{Sl))r+il2{{SlS2))r 

+6(A;2 - k^){{SiSl))r + 3(Zo - 2h + h){{SlSl))r 

+ h{{SiS2S^))r + 3(^4 - l2){{SlS2Sl))r 
+ 3{ki-2k2 + h){{SiSiSl))r 

+ {-lo + 3/i - 3^3 + k){{SlSlSi))r, (11) 

We note that, for the sake of simplicity, the superscript z is omitted from the correlation functions on the right-hand 
side of Eq. (Ill|) . The coefficients in Eq. (fTTj) are defined as follows: 

h = cosh(JV)F(x)U=o, ki - sinh(JV)F(a;)|,^o, 

h = sinh2(JV)F(x)U^o, h = cosh(JV)sinh(JV)F(a;)U=o, 

^3 = cosh2(JV)F(a;)U=o, ^3 = smh^{jy)Fix)U=o, 

h = cosh(JV)sinh^(JV)i^(a;) 1:^^=0, k^ = cosh^(JV)sinh(JV)F(a;)|2;=o, 

h = cos\i\jV)F{x%=o. (12) 



Next, the average value of the perimeter spin in the system can be written as follows, and it is found as 

™i = {{Ss))r = ((1 + 5'(5sinh(JV) + (^o')'{cosh(JV) - l}))rF{x + 7)U=o, (13) 



with the coefficients 



{{S^))r = ai (1 - {{{S^r))r) + a,{{S^))r + aMS'or))r, 



ai = F(7), 

02 = sinh(JV)F(a; + 7)|:c=o, 

as = cosh(JV)F(x + 7)|x=o, 



(14) 



(15) 



where 7 = ((? — 1)^ is the effective field produced by the (q — 1) spins outside the system and A is an unknown 
parameter to be determined self-consistently. In the effective-field approximation, the number of independent spin 
variables describes the considered system. This number is given by the relation v = (((-Sf )^'^))r- As an example for 
the spin-1 system, 2S = 2 which means that we have to introduce the additional parameters (((»S'o)^))r and (((S'|)^))r 
resulting from the usage of the Van der Waerden identity for the spin-1 Ising system. With the help of Eq. ([5]), 
quadrupolar moment of the central spin can be obtained as follows 



(((^o)')). = ( ( n [1 + ^|sinh( JV) + (5|)2{cosh(JV) - 1}] 



Gix)\ 



where the function G{x) is defined as 



Gix)= I dD,P{Di)g{x,Di). 



(16) 



(17) 



Definition of the function g{x,Di) in Eq. (IT71) is given as follows and the expression in Eq. ^TE\\ can be evaluated by 
using the eigenvalues and corresponding eigenvectors of the effective Hamiltonian matrix in Eq. ([7]) . 



gixjDi) 



E„=iexp(/3A„)„^^ 
2 cosh(/3x 



^(^„|(5f)^|^„)exp(/3A„) 



(18) 



2cosh(/3x) +e-f^^-' 
Hence, we get the quadrupolar moment by expanding the right-hand side of Eq. (|16p 

(((^o)'))r = ro + 3ni{{Si))r + 3{ri~ro){{Sl))r + 3r2{{SiS2))r 

+6(n2 - n,){{SiSi))r + 3(ro - 2n + r3){{SfSl))r + n3((^i^2^3)). 
+3(r4 - r2){(SiS2Sl))r + 3(ni - 2^2 + n4)((5i5|^l))r 



+ (-ro + 3ri - 3r3 + r5){{SlSlSl)) 



■il iri 



(19) 



with 



ro = G(0), 

Ty =cosh(JV)G(a;)|^=o, 

r2=sinh2(JV)G(x)U=o, 

rz = coaY? {JV)G{x)\x=o, 

Ti = cosh(JV)sinh^(JV)G(a;)U=o, 

r5 = cosh'(JV)G(:E)U=o, 



m =sinh(JV)G(a;)|2.^o, 

n2 = cosh(JV)sinh(JV)G(x)|j;=Oj 

^3 = sinh^(JV)G(a:)|a;=o, 

Ui = cosh^(JV)sinh(JV)G(x)U=o, 



Corresponding to Eq. ((13)) . 

mif))r = ((1 + ^o'sinh(JV) + (5o^)2{cosh(JV) - l}))rG{x + 7), 



(20) 



(21) 



{{Sl))r ^b,{l~ {{{S^f))r) + b2{{S^,))r + ^{{{3^)')^ 



(22) 



where 

bi = G(7), 

62 = sinh(JV)G(x + 7)U=o, 

63 - cosh(JV))G(a; + 7)U=o. (23) 

Eqs. (|lip . ([T4)) . (fT9|) and (l22t are the fundamental spm identities of the system. When the right-hand sides of Eqs. 
((HI) and ()16|) are expanded, the multispin correlation functions appear. The simplest approximation, and one of the 
most frequently adopted is to decouple these correlations according to 

{{snS^r...S^))^ - ((5f )), {{{S^r))^ ... ((^f )). , (24) 

for i y^ j ^ ... 7^ &. The main difference of the method used in this study from the other approximations in 
the literature emerges in comparison with any decoupling approximation (DA) when expanding the right-hand sides 
of Eqs. (|8]) and (|16p . In other words, one advantage of the approximation method used in this study is that no 
uncontrolled decoupling procedure is used for the higher-order correlation functions. 

For spin-1 Ising system with q = 3, taking Eqs. pT|) . ((T4)) . p9)) and (|22|) as a basis, we derive a set of linear equations 
of the spin identities. At this point, we assume that (i) the correlations depend only on the distance between the 
spins, (ii) the average values of a central spin and its nearest-neighbor spin (it is labeled as the perimeter spin) are 
equal to each other with the fact that, in the matrix representations of spin operator S, the spin-1 system has the 
properties (Sg)^ — 5| and (5*1)^ = (5*1)^. Thus, the number of the set of linear equations obtained for the spin-1 
Ising system with q = 3 reduces to twenty one, and the complete set is given in Appendix IB] 

If Eq. (jBip is written in the form of a 21 x 21 matrix and solved in terms of the variables Xi[{i = 1, 2, ..., 21){e.g., xi — 
(('S'o))r: 2:2 = ((S'iS'o))r, ■•■)] of the linear equations, all of the spin correlation functions, as well as magnetizations and 
quadrupolar moments can be easily determined as functions of the temperature and Hamiltonian parameters. Since 
the thermal and configurational average of the central spin is equal to that of its nearest-neighbor spins within the 
present method then the unknown parameter A can be numerically determined by the relation 

{{S^))r ^ {{Si))r or a;i=X4. (25) 

By solving Eq. (j25p numerically at a given fixed set of Hamiltonian parameters we obtain the parameter A. Then 
we use the numerical values of A to obtain the spin correlation functions which can be found from Eq. (|B1|) . Note 
that A = is always the root of Eq. (pS)) corresponding to the disordered state of the system. The nonzero root of A 
in Eq. (1^ corresponds to the long-range ordered state of the system. Once the spin identities have been evaluated 
then we can give the numerical results for the thermal and magnetic properties of the system. Since the effective 
field 7 is very small in the vicinity of ksTc/J, we can obtain the critical temperature for the fixed set of Hamiltonian 
parameters by solving Eq. ([25]) in the limit of 7 ^' then we can construct the whole phase diagrams of the system. 
Depending on the values of Hamiltonian and crystal field distribution parameters, there may be two solutions [i.e., 
two critical temperature values which satisfy Eq. (1251) ] corresponding to the first (or second) and second-order phase- 
transition points, respectively. We determine the type of the transition by looking at the temperature dependence of 
magnetization for selected values of system parameters. 

III. RESULTS AND DISCUSSION 

In this section, we will discuss the effect of the crystal field distributions defined in Eqs. © and (|3]) on the global 
phase diagrams of the system where the second and first order transitions are shown by solid and dashed curves, 
respectively with tricritical points (shown by hollow circles) for honeycomb (q = 3) and square (q = 4) lattices. Also, 
in order to clarify the type of the transitions in the system, we will give the temperature dependence of the order 
parameter. 

A. Phase diagrams of the system with dilute crystal field 

In this section, we illustrate the phase diagrams and magnetization curves of the system with a dilute crystal field 
distribution defined in Eq.(l2]) where crystal field D is turned on, or turned off with probabilities p and (1 — p) on the 
lattice sites, respectively. In Figs. ([T^) and (Ht), we plot the phase diagrams of the system in {ksTc/J — D/J) plane 
for honeycomb and square lattices with coordination numbers q — 3 and g = 4, respectively. As seen in Figs. ((1^) 



and ^p), phase diagrams of the system can be divided into three parts with different concentration values p. For 
the curves in the first group with p < p*, the system always exhibits a second order phase transition with a finite 
critical temperature ksTc/ J which extent to D/ J — s- — oo. If the concentration p reaches its critical value p* then the 
critical temperature depresses to zero. Physical reason underlying this behavior can be explained as follows; when we 
select sufficiently large negative crystal field values {i.e.D/J -^ — cxd), all of the spins in the system tend to align in 
5 = state. As p increases starting from zero, the ratio of spins which aligned in 5 = state increases, and therefore, 
magnetization weakens, and accordingly, critical temperature of the system decreases. According to our numerical 
results, the critical concentration value is obtained as p* = 0.3795 for q = 3 and p* = 0.5875 for q — A. In the second 
group of the phase diagrams in Figs. (H^) and ^), the system exhibits a reentrant behavior of second order, whereas 
the curves in the third group, exhibit a reentrant behavior of first order with a tricritical point at which a first order 
transition line turns into a second order transition line. Besides, the curves which exhibit a reentrant behavior of 
first (or second) order, depress to zero at three successive values of crystal field D/J = —3.0, —2.0, —1.0. Moreover, 
in D/J -^ oo limit, the system behaves like spin-1/2 for p — 1.0. In the case of p 7^ 0, the ratio of spins that behave 
like S = ±1 increases as p increases. Therefore, for < p < 1.0, all transition lines have finite critical temperatures 
which increase with increasing p values for D/J -> cxa. At this point, we also note that if we select D/J = in Eq. 
([2]), all lattice sites expose to a crystal field Di/J = independent from p. Hence, all transition lines intersect each 
other on the point {D/J^ksTc/J) = (0,1.3022) for <? = 3, and {D/J^ksTc/J) = (0,1.9643) for q = 4. Meanwhile, 
previous studies based on EFT are not capable of obtaining first order transition lines of the system. From this point 
of view, we see that our method improves the results of the other EFT works and we take the conventional EFT 
method one step forward by investigating the global phase diagrams, especially the first-order transition lines that 
include reentrant phase transition regions. 

TABLE I: Critical concentration p* obtained by several methods 
and the present work for honeycomb (g = 3) and square [q = 4) 
lattices. 



Lattice EFT-ia>2 EFT-HiS MFT^^^ PA^S Present Work 
g = 3 0.484 0.492 LO 0.5 0.3795 

q = 4 0.604 0.610 1.0 0.667 0.5875 



On the other hand, Figs. ([T]d) and ([TJl) shows the phase boundary in [ksTc/ J ~ p) plane which separates the 
ferromagnetic and paramagnetic phases with D/J — > — cxd. According to this figure, critical temperature ksTc/J of 
system decreases gradually, and ferromagnetic region gets narrower as p increases, and kgT^/J value depresses to 
zero at p = p* . Such a behavior is an expected fact in dilution problems. Numerical value of critical concentration 
p* for honeycomb {q = 3) and square {q = 4) lattices is given in Table HJ and compared with the other works in the 
literature. As seen in Table HI numerical values oi Pc for q = 3 and q = A are new results in literature. Furthermore, 
MET— kiS- predicts that the system always has a finite critical temperature and exists in a ferromagnetic state at lower 
temperatures in D/J -^ —00 limit, except that p — 1.0. This artificial result can be regarded as a failure of the MET. 

In Fig. ([2]), we plot the temperature dependencies of magnetization curves corresponding to the phase diagrams 
depicted in Fig. ([T]) for q = 3. As seen in Fig. ([2]), as p increases then critical temperature ksTc/J values decrease for 
D/J < 0, except the reentrant phase transition temperatures which occur at low temperatures. On the other hand, 
effect of increasing p values on the shape of magnetization curves depends on value of D/J. Namely, in Figs. ([2^) and 
(HJd) we see that ground state saturation values of magnetization curves decreases as p increases for D/J = —10.0 and 
—3.1. Moreover, for D/J — —3.1, magnetization curves of the system exhibit a broad maximum at low temperatures 
for p = 0.37, and a reentrant behavior of second order for p = 0.4. If we select D/J = —2.5 as in Fig. ([It), saturation 
values of magnetization curves remain unchanged for p — 0,0.2,0.3 and tend to decrease for p > 0.3. li p increases 
further, a reentrant behavior of second order appears for p = 0.53, and we see a broad maximum at low temperatures 
for p = 0.517 and 0.5 which tends to depress as p decreases. This broad maximum behavior of magnetization curves 
originates from the increase in the number of spins directed parallel to the z-direction with increasing temperature, due 
to the thermal agitation. For D/J = —2.0 in Fig. (IlJl), magnetization curves saturates at m = 1 at the ground state 
and reentrant behavior disappears. If we increase D/J further, for example for D/J = —1.5 (Fig. (dh)), another type 
of reentrant behavior occurs in the system in which a first order transition is followed by a second order transition. 
Finally, for sufficiently large positive values of crystal field, magnetization curves always saturate at to = 1 and the 
system always undergoes a second order phase transition from a ferromagnetic phase to a paramagnetic phase with 
increasing temperature, which can be seen in Fig. ([2f) with D/J = 10.0. As a common property of the curves in 
Fig. ([5]), we see that effect of p on the saturation values, as well as temperature dependence of magnetization curves 
strictly depend on the strength of D/J. Hence, according to us, the presence of dilute crystal fields on the system 



should produce a competition effect on the phase diagrams of the system. We also note that, although it has not been 
shown in the present work, magnetization curves for g = 4 corresponding to the phase diagrams depicted in Fig. ([TJ;) 
exhibit qualitatively similar behavior with those of Fig. ([2]) with q — 3. 

As seen in Fig. ([T]), for a dilute crystal field distribution defined in Eq. ^, the global phase diagrams which are 
plotted in (kBTd J — D / J) plane, as well as the phase boundaries in {kBTc/ J — p) plane for q — Z exhibit qualitatively 
similar characteristics when compared with those for (7 = 4. Hence, in order to examine the phase diagrams which are 
plotted in {kBTc/ J — D/ J) plane in Figs. ([Iji) and ([TJ;) in detail, we plot the evolution of the global phase diagrams 
in Fig. ([31) only for q — A. From this point of view, Fig. (|3^) shows how the phase diagrams in Fig. ([T};) evolve 
when the concentration p changes from 0.5 to 0.6. As seen in Fig. (|3^), we observe a second order phase transition 
line with a finite critical temperature /cbTc/ J which extent to D/ J — > — oo for p — 0.575. If p increases, namely 
for p = 0.580 and 0.583, we see that a low temperature transition line arises between —4.0 < D/J < —3.0, as well 
as a high temperature phase boundary which extents to D/J -^ — cx). If p increases further, such as for p = 0.584, 
0.585 and 0.587, high temperature phase boundary is gradually connected to the transition line which arises between 
—4.0 < D/J < —3.0, and the phase diagrams exhibit a bulge on the right hand side of {kBTc/ J — D/J) plane, whereas 
another transition line emerges within the range of — oo < D/J < —4.0, which disappears for p > 0.587. Similarly, 
evolution of the phase diagrams in Fig. ([!];) when the concentration p changes from 0.6 to 0.7 can be seen in Fig. 
([3)d). As seen in this figure, the curves for p = 0.60, 0.62, 0.64, 0.66 exhibit a reentrant behavior of second order, 
while for p — 0.68 reentrance disappears and for p ~ 0.70 and 0.71 double reentrance with three successive second 
order phase transitions occurs in a very narrow region of D/J. On the other hand, increasing values of p generates 
first order phase transitions with tricritical points, as well as reentrant behavior of first order. This phenomena is 
illustrated in Fig. ^jp). From Fig. (jSt), we see that, the second order transition temperatures decrease as absolute 
value of D/J increases, and turn into first order transition lines at tricritical points. Evidently, the phase diagrams 
change abruptly for p > 0.7164. Hence, the behavior of the p = 0.7163 curve is completely different from that of 
p = 0.7164. Namely, first order transition temperatures of the system for p < 0.7163 and p > 0.7164 depress to zero 
at D/J = —3.0 and D/J = —2.0, respectively. In order to investigate the phase transition features of the system 
further, we should continue increasing the value of p. In Fig. ^JX), we see that the curves for p = 0.72, 0.74, 0.76 
and 0.78 exhibit a reentrant behavior of first order, whereas the curves with p — 0.80, 0.82, and 0.84 exhibit double 
reentrance with two first order and a second order transition temperature. 

B. Phase diagrams of the system with random crystal field 

Next, in order to investigate the effect of the random crystal fields defined in Eq.Q on the thermal and magnetic 
properties of the system, we represent the phase diagrams and corresponding magnetization curves for honeycomb 
{q = 3) and square lattices {q = 4) throughout Figs. ^ and ([S]). 

We note that random crystal field distribution given in Eq. ([3]) with D/J = corresponds to a bimodal distribution 
function, while for A/ J = 0, we obtain pure BC model with homogenous crystal field D/J. In Fig. ^, phase diagrams 
of the system corresponding to the bimodal distribution function are shown in {kBTc/ J — A/ J) plane. For a bimodal 
distribution, the phase diagrams have symmetric shape with respect to A/ J which comes from the fact that p = 1/2, 
and as seen in Fig. ([4]), transition temperatures are second order, and it is clear that the system exhibit different 
characteristic features depending on the coordination number q. Namely, for 9 = 3, transition temperature decreases 
with increasing A/J and exhibits double reentrance with three second order phase transition temperatures, then 
falls to zero at A/J = 3.0 (left panel in Fig. ^). On the other hand, as seen on the right panel in Fig. ^, as 
A/J increases then the transition temperature of the system for g = 4 decreases and remains at a finite value for 
A/J — > oo which means that ferromagnetic exchange interactions ior q — 3 are insufficient for the system to keep its 
ferromagnetic order for A/J > 3.0, while for q = A these interactions are dominant in the system, and the presence 
of a disorder in the crystal fields cannot destruct the ferromagnetic order. 

At the same time, in order to see the effect of the random crystal fields with A/J, D/J ^ on the phase diagrams 
and magnetization curves of the system for q — 3 and 4, we plot the phase diagrams in {kBTc/ J — D/J) plane in 
Fig. ([5]) and variation of the corresponding magnetization curves with temperature in Figs. © and ([7]), respectively. 
At first sight, it is obvious that the phase diagrams in Fig.® represent evident differences in qualitative manner 
with coordination number q. That is, as seen in Fig. ([5^), the curve corresponding to A/J = represents the 
phase diagram of pure BC model for a honeycomb lattice which exhibits a reentrant behavior of first order with first 
and second order transition lines, as well as a tricritical point. From Fig. ([5^), we see that as A/J increases then 
the tricritical point and first order transitions disappear, and the first order reentrance turns into double reentrance 
with three transition temperatures of second order, and phase transition lines shift to positive crystal field direction 
without changing their shapes. On the other hand, the situation is very different for a square lattice. Namely, as 
seen in Figs, ^li.) and (tSb); A/J = curves for g = 3 and g = 4 are qualitatively identical to each other. However, 
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as seen in Fig. ([5]d), for A/ J ^ 0, first order transition lines and tricritical points do not disappear from the system 
for g = 4, but shift to negative crystal field values. Besides, the system does not exhibit double reentrance for q — 4:. 
Furthermore, for A/J > 1.5 in Fig. ([5^), and A/J > in Fig. ([Sb), all phase diagrams exhibit similar behavior as 
D/J varies. Namely, critical temperature ksTc/J in Fig. ([5^) reduces to zero at D/J = A/J — 3.0. On the other 
hand, first order transition temperatures in Fig. ^}p), reduce to zero at D/J = —A/J. 

It is important to note that these observations are consistent with the results shown in Figs. ((T^) and ((TJ;). In 
other words, the distribution function given in Eq. ^ with p = 0.5 and D/J = 2Do/J is identical to Eq. Q for 
A/J = ±Dq/J and D/J — Dq/J. For example, according to Eq. ^, if we select Dq/J = 4.0 with p = 0.5, it means 
that half of the spins on the lattice sites expose to a crystal field D/J — 0, while a crystal field given by D/J = 8.0 
acts on the other half of the spins. On the other hand, if we select A/J = ±Dq/J and D/J — Dq/J by using Eq. (O, 
we generate the same distribution again. Hence, we expect to get the same results in Figs. ([T|) and (0 for these system 
parameters. For instance, for Dq/J — 4.0 in Fig. (H^), we get D/J — 8.0, and the system exhibits a ferromagnetic 
order in the ground state, which can also be seen in Fig. ([5^) with a critical temperature ksTc/J — 1.4395. These 
conditions are also valid for q — A, and for the whole temperature region on the phase diagrams. Therefore, the state 
(para-or ferro), as well as thermal and magnetic properties of a selected {ksT/ J,D / J) point with respect to p — 0.5 
curves in Figs. ([T^) and ([T];) is identical to the state of a point (ksT/J, D/2J) in Figs. ([5^) and ([SJd) with respect to 
the curve A/J = ±D/J, respectively. Moreover, the qualitative differences between Figs. ([5^) and ([SId) mentioned 
above are strongly related to the percolation threshold value of the lattice. Namely, distribution function Eq.Q is 
valid only for p — 0.5. However, as seen in Table HI we obtain pc < 0.5 for g = 3, and Pc > 0.5 for q = 4. 

In Fig. (|n]), we examine the temperature dependence of magnetization curves for g = 3, corresponding to the phase 
diagrams shown in Fig. ([S^) with D/J = —1.0. Fig. (H^), shows how the temperature dependence of magnetization 
curves evolve when A/J changes. According to Fig. ([B^), magnetization curves saturate at a partially ordered state 
at low temperatures. Besides, for A/J = 1.68, 1.74 and 1.80, the system undergoes three successive phase transitions 
of second order, which confirms the existence of double reentrance. Similarly, Fig. (JBJd) shows how the shape of the 
magnetization curves change as D/J changes for constant A/J ~ 6.0. As seen in Fig. (Hb), magnetization curves 
exhibit a second order phase transition from a ferromagnetic (fully ordered) to a paramagnetic phase at certain values 
of crystal field, namely at D/J — 4.0, 4.1 and 4.4, whereas for D/J — 3.3, 3.5, 3.7 and 3.9 the system can only achieve 
a partially ordered phase. In addition, the curves corresponding to D/J ~ 3.5, 3.7 and 3.9 exhibit a broad maximum 
at low temperatures, and then decrease as the temperature increases, whereas for D/J = 3.3, we observe double 
reentrance. Fig. ([7]) shows the magnetization curves for g = 4, corresponding to the phase diagrams shown in Fig. 
(I5}d). In Fig. d?^), we see that magnetization curves exhibit a second order phase transition from a paramagnetic 
phase to a fully ordered ferromagnetic phase for A/J = and 5.0. On the other hand, the curves corresponding to 
D/J = 5.5, 5.8 and 6.0, saturate at a partially ordered state at low temperatures, and exhibit a broad maximum with 
increasing temperature which depresses gradually as A/ J increases, then fall rapidly at a second order phase transition 
temperature. The broad maximum behavior observed in these curves disappears for A/J = 8.0. Additionally, Fig. 
([7)3) represents the magnetization versus temperature curves for g = 4 with with D/J — —4.0. In Fig. ((Tb), it 
is clearly evident that, at low temperatures, the system saturates at a partially ordered phase for A/J = 4.0, 5.0, 
6.0 and 7.0, while for A/J = 3.7 and 3.8, a reentrant behavior of first order occurs. Again we see that, there is a 
competition between ferromagnetic exchange interactions and disorder effects in crystal fields which determines the 
saturation values and temperature dependence of magnetization curves of the system. 

Finally, dependence of magnetization of the system on the crystal field A/J for fixed temperature values ksT / J ~ 
0.01,0.05,0.1 and 0.2 with D/J = 2.0 is shown in Fig. ([5]) for g = 3 and 4, respectively. We see that at sufficiently 
low temperatures such as ksT/J = 0.01, magnetization curves exhibit three phases for g = 3. A first order transition 
is characterized by a gap in this figure. On the left panel in Fig. ([S]), which is plotted for g = 3, we observe two 
successive first order phase transitions for ksT/J = 0.01. The first one is from the fully ordered ferromagnetic phase 
{m = 1.0) to the partly ordered phase (to = 0.47), and the other is from partly ordered phase to disordered phase 
(to = 0.0). On the other hand, according to the right panel in Fig. ([8]), the system can not reach a paramagnetic 
phase at the ground state for g = 4. Hence, we observe two phases. Namely, for ksT/J = 0.01, a first order phase 
transition from a fully ordered phase (m = 1.0) to a partly ordered phase (m = 0.59). Then, saturation magnetization 
of the partly ordered phase reduces continuously to (to, — 0.527). Moreover, the first order transitions disappear with 
increasing temperatures, both for g = 3 and 4. Since in Figs. ([5^) and (03), all phase diagrams exhibit similar 
behavior as D/J varies, behavior of magnetization curves in Fig. ([S]) should be the same for different D/J values. 
Namely, the left panel of Fig. ([5]) can be regarded as the variation of magnetization with A/J within the range 
D/J + 1.0 < A/J < D/J + 4.0 for g = 3. for a selected value of D/J. It is possible to observe the similar behavior 
on the right panel in Fig. (jH)). This general behavior is an expected result, since the behavior of the system depends 
on the relation between A/J and D/J parameters for a given temperature, not their values independently. 



IV. CONCLUSIONS 

In this work, we have studied the phase diagrams of a spin-1 Blunie-Capel model with diluted and random crystal 
field interactions on two dimensional lattices. We have introduced an effective-field approximation that takes into 
account the correlations between different spins in the cluster of a considered lattice and examined the phase diagrams 
as well as magnetization curves of the system for different types of crystal field distributions, namely, dilute crystal 
fields and a double peaked delta distribution, given by Eqs. (I2|) and ([3]), respectively. 

For dilute crystal fields, we have given a detailed exploration of the global phase diagrams of the system in ksTc/ J — 
D/ J plane with the second and first order transitions, as well as tricritical points. We have also shown that the system 
with dilute crystal fields exhibits a percolation threshold value Pc which can not be predicted by standard MFA. In 
addition, we have observed multi-reentrant phase transitions for specific set of system parameters. 

On the other hand, we have investigated the effect of the random crystal field distribution characterized by two 
crystal field parameters D / J and A/J on the phase diagrams of the system. As a limited case, we have also 
focused on a bimodal distribution with D/ J — 0. Particulary, we have reported the following observations for a 
bimodal distribution: It has been found that the phase diagrams have symmetric shape with respect to A/J which 
comes from the fact that p = 1/2. The transition temperatures are of second order, and the system exhibit different 
characteristic features depending on the coordination number q. Besides, we have realized that the system may exhibit 
clear distinctions in qualitative manner with coordination number q for random crystal fields with A/J, D / J ^ 0. 
Moreover, we have discussed a competition effect which arises from the presence of dilution, as well as random crystal 
fields, and we have observed that saturation values of the magnetization curves are strongly related to these effects. 

As a result, we can conclude that all of the points mentioned above show that our method improves the conventional 
EFT methods based on decoupling approximation. Therefore, we hope that the results obtained in this work may be 
beneficial from both theoretical and experimental points of view. 
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Appendix A: Fundamental correlation functions of the system for a square lattice 



Magnetization of the central spin for a square lattice is given as follows 



((^o)) 



Mo + 4ci((S'i))^ + 4(/i2 - ^^o){{Sf))r 

+(i^^l{{SlS2))r + 12(C2 - C,){{SiSi))r 

+6(^10 - 2^2 + M3)((^?5|))r + 4c3((5l52^3))r 

+ 12(^^4 - fll){{SiS2Sl))r + 12(C4 - 2C2 + Ci ) ((^i 5|^|)). 

+4(^5 - 3^3 + 3^2 - Aio)((<5'l'S'|5'|)).r + ^J.s{{SlS2S3S4))r 

+4(C5 - C3){{SiS2S3Sl))r + Q{^il ~ 2/^4 + Me) ((5l^2^3'^4 ))r 

+4(C6 - 3C4 + 3C2 - Ci){{SiSlSlSl))r 

+ (/iO - 4^2 + 6/i3 - 4/i5 + li7){{SlSlSlSl))r. 



where the coefficients are given by 



Mo = 


= FiO), 




A^i = 


= sinh2(JV)F(x)U=o, 


Cl 


m '- 


= cosh(JV)F(x)U=o, 


C2 


Ma = 


= cosh2(JV)F(x)U=o, 


C3 


M4 = 


= sinh2(JV)cosh(JV)F(a;)U=o, 


C4 


Ai5 = 


= cosh3(JV)F(x)U=o, 


C5 


Me = 


= sinh2(JV)cosh2(JV)F(x):,=o, 


ce 


M7 = 


= cosh4(JV)F(x)U^o, 




M8 = 


= sinh4(JV)F(a;),=o. 





sinh(JV)F(x):,=o, 
sinh(JV) cosh(JV)F(x)2;=o, 
sinh3(JV)i^(x),=o, 
cosh^(JV) sinh( JV)F(x)a;=o, 
sinh^( JV) cosh(JV)F(a:)x=o, 
cosh^(JV) sinh( JV)F(a::)^^o, 



Quadrupolar moment corresponding to equation (jBip defined as 

(((^o)')) = Po + 4r;i((5i)),+4(M2-Mo)((5?)). 

+ 6pi{{SlS2))r + 12irj2 - 7n){{SlSi))r 

+ 6(mo - 2p2 + P3){{S!Sl))r + 4773((^1^2^3)). 

+12(m4 - Mi)((5i^2^|))r- + 12(^4 - 2772 + Vi){{SiSiSi))r 

+4(m5 - 3p3 + 3p2 - P0){{SfSlSl))r + P8{{SlS2S3Si))r 
+4fe - ^3)((5l52^3^4»r + 6(pi - 2m4 + Pe) ((^l ^2^3 ^4 )) 
+4(776 - 37^4 + 3772 - T^i){{SiSlSlSl))r 

+ (P0 - 4/92 + 6p3 - 4M5 + P7){{SfS^SiSl))r, 



where 



= G(0), 

:sinh2(JV)G(a;)U=o, 
= cosh(JV)G(a;)U=o, 
cosh2(JV)G(a;) 



\x=0, 



Po 

Pi 

P2 
P3 ^ 

M4 = sinh^(JV)cosh(JV)G(a;)U=o, 

M5=cosh3(JV)G(a;)U=o, 

pe ^ sinh^(JV)cosh^(JV)G(x)^=o, 

M7 = cosh4(JV)G(a;)U=o, 

Ms = sinh''(JV)G(x)j;=o- 



sinh(JV)G(a;)^=o, 

sinh( JV) cosh( JV)G(x)a;= 

sinh3(JV)G(a;) 



)x=0, 



m 
m 
m 

rii = cosh^(JV)sinh(JV)G(a;), 

775 = sinh''^(JV)cosh(JV)G(a:;), 

776 = cosh^(JV)sinh(JV)G(a:;), 



=0, 
=0, 
=0, 



Finally, perimeter spin identities are as follows 

a,{l-{{{S,f))r 

UJl +W2((S'o))r + 






)+«2((^o)),.+«3(((5o)'))., 
(W3 ~UJi){{Sl))r 



(Al) 



(A2) 



(A3) 
(A4) 
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with the coefBcients 

ai = F(7) LOi = G(-f)\x=o 

a-2. = sinh( JV)F(x + 7) ^2 = sinh(JV)G(a; + 7)^=0 

as — cosh( JV)i^(x + 7) cja = cosh( JV)G'(x + 7)|a:=o 

where 7 = (g — 1)A with g = 4, and the functions F[x) and G(x) are defined in equations ([9]) and (|T7)) . 
Appendix B: The complete set of twenty one linear equations of a honeycomb lattice 



+ 6(fe - kr)l,l,SxSl))r + 3(^0 - 2/1 + h){.{.SlSl))r 
+ h{{S,S2S3))r + m - l2){{SlS2Sl))r 

+3iki-2k2 + h){{SiSiSl))r 

+{-io + ill - 313 + i5){{s!sisl))r 

{{SlSo))r = (3^1 - 2lo){{Si))r + 3ki{{Sf))r + Hk " 2^1 + ^2 + ^3)((^l5|)). 
+ 6(^2 - ki){{SlSi))r + k3{{SiS2Sl))r 

+{-lo + 'ill - ih - 3h + ih + h){{SiSlSl))r 

+ 'i{ki-2k2 + ki){{SlSlSl))r 
{{SlS2So))r = ilo-3h+3l2+3l3){{SiS2))r + {Gk2-3ki){{SiSl))r 
+ i-lo + 3/i - 3^2 - 3^3 + 3^4 + h){{SlS2Sl))r 

+ {3ki - 6/c2 + fcs + 3k4){{SiSiSl))r 
{{Si))r = ai{l~{{{S^f))r)+a2{{S^))r+a,{{iS^)')),. 
{{SiS2))r = ai((^i)),. + a2((5'o^i». + (a3-ai)((^i^o))r- 
{{SiS2S3))r = ai((S'iS'2))r + a2(('S'o5'iS'2))r + (03 - ai)((<5'iS'2S'o))r 

{{Sl))r = biil~{{iS^f))^)+b2{{S^))r+h{{iS^r))r 

{{SiSl))r = bi{{Si))r+b2{{SQSi))r + {b3-bi){{SiS^,))r 

{{SfSl))r. = bi{{Sl))r+b2{{SoSl))r + ib3-bi){{SlS^o))r 

{{SoSf))r = b3{{So))r+b2{{S^))r 

((SoSiSi))r = b3{{SoSi))r+b2{{SiS^))r 

{{SoSfSl))r - b3{{SoSf))r+b2{{SfS^o))r 

{{SiS2Sl))r = bi{{SiS2))r+b2{{SoSiS2))r + ib3-bl){{SiS2S^))r 

((SiSiSl))r = bi{{SiSl))r + b2{{SoSiSl))r + (63 - bi){{SiSlS^))r 

{{SfSiSl))r = bi{{SfSi))r + b2{{SoSfSi))r + (63 - bl){{SfSlS^o))r 

{{S^))r = ro + 3ni((^i)), + 3(ri-ro)((^?)). + 3r2((5i^2))r 

+6(n2 - ni){(SiSl))r + 3(ro - 2ri + r3){{SfSl))r + n3((5i52^3))r 
+3(r4 - r2)((5i52^|))r- + 3(ni - 2n2 + ni){{SiSiSi))r 
+ (-ro + 3ri - 3r3 + r5){{SfSlSl))r 
{{SiS^))r = (3ri - 2ro)((5i)), + 3ni{{Sl))r + (3r2 + 3ro - Gn + 3r3){{SiSl))r 
+6(^2 - ni){{SfSl))r + n3{{SiS2Sl))r 
+ (-ro + 3ri - 3r2 - Sr^ + Sn + r5){{SiSiSl))r 
+3(ni - 2712 + n4){{SfSiSl))r 
{{SfS^))r = (3ri - 2ro)((52)),, + 37ii((5i)), + (3r2 + 3ro - Gn + 3r3)((5252'))r 
+6(n2 - ni)((S'iS'2)),- + (3ni - 6n2 + ^3 + 3714) ((S"! 6*2 S'f))r 
+ (-ro + 3ri - 3r2 - 3r3 + 3r4 + r5){{SfSiS^3))r 
{{SiS2S^))r = iro-Sri+3r2 + -ir3){{SiS2))r + i-Sni + 6n2){{SiSl))r 



12 

(-ro + 3ri - 3r2 - Sra + 3r4 + r5){{SiS2S'^))r 

+ (3ni - 6n2 + ns + 3n4)((5iS'2S'f ))r 
((^i52^o'))r = (ro-3ri+3r2 + 3r3)((5i^2))r + (-3ni+6n2)((5i^2))r 

(-ro + 3ri - 3r2 - 'ir^ + Sr^ + r5)((S'iS'|S'|))r 

+ (3ni - 6n2 + ns + 3n4)((S'i5'25'|))r 
(('5i^|>5o))r = (ro - 3ri + 3r2 + 3r3)((525|)), + (-3ni + 6n2){{SiSl))r 

+ (3ni - 6n2 + ns + 3n4)((5iS'2S'|))r 

+ (-ro + 3ri - 3r2 - 3r3 + 3r4 + r5)((^?52'^3'))r- 

(Bl) 
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FIG. 1: (a) Phase diagrams of the system for g = 3 in a {kBTc/J — D/J) plane corresponding to dilute crystal field distribution 
defined in Eq. ([2]). The solid and dashed lines correspond to second- and first-order phase transitions, respectively. The open 
circles denote the tricritical points, and the numbers on each curve represent the value of concentration p. (b) Phase diagrams 
of the system for g = 3 in a {kBTc/J — p) plane with a selected value of the crystal field D/J = —10.0. (c) Phase diagrams of 
the system for g = 4 in a {ksTc/J — D/J) plane corresponding to dilute crystal field distribution defined in Eq. ((2|. The solid 
and dashed lines correspond to second- and first-order phase transition, respectively. The open circles refer to the tricritical 
points, and the numbers on each curve represent the value of concentration p. (A) Phase diagrams of the system for g = 4 in a 
(kBTc/J — p) plane with a selected value of the crystal field D/J — —10.0. 
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FIG. 2: (Color online) Temperature dependence of magnetization corresponding to Fig. ([TJi) with some selected values of 
crystal field, (a) D/J = -10.0, (b) D/J = -3.1, (c) D/J = -2.5, (d) D/J = -2.0, (e) D/J = -1.5, and (f) D/J = 10.0. The 
numbers on each curve denote the value of concentration p. The solid and dashed lines correspond to second- and first-order 
phase transitions, respectively. 



16 












^ 


■ 






^^^ 


^=^ 


0.66 


^ 


^F^ 




0.64 


\yZ 


t 


/ 




■ 0.62 \ 

;o.6o \ y 


W/ 


^p=0.71 




■ \}v 


/// 


^ 


^0.7 




■ /// 


(L 


^ 


^0.68 




■ ^4 


\M 






q=4 



D/J 



D/J 



(c) 




0.7164 



(d)ia 

1.4 

1.2 

1.0 
-J 
O 0.8 

CO 
^ 0.6 

0.4 

0.2 

0.0 



0.78 ^....=:^ 


^^^ 


076 \ ^.^:s^^^^^ 




74 \ \s^^^^^^^ 




"p=0.72 \^>^^^§C<i 




/wvv/kT ^-^^ 




- //Ilfif^^ °-^^ 




■ I I ^ \\'^ ' 0.8 

\ S N ^ ^ ^ ^ ^ 

"" ~ - - ■- ^ __^ \ y 

1 . 1 . ~ T - -T---i^i- .1 


q=4 



-2.0 



-2.6 



-1.8 



-1.6 



D/J 



D/J 



FIG. 3: (Color online) Evolution of the phase diagrams corresponding to Fig. Ht). The numbers on each curve denote the 
value of concentration p. The solid and dashed lines correspond to second- and first-order phase transitions, respectively. The 
open circles indicate the tricritical points. 
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FIG. 4: Phase diagrams of the system in a (fesTc/ J — A/J) plane for a bimodal crystal field distribution corresponding to Eq. 
(|3]) with D I J — 0.0. Left and right-hand side panels are plotted for g = 3 and g = 4, respectively. 
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FIG. 5: (Color online) Phase diagrams of the system in a {ksTcj 3 — D/J) plane corresponding to random crystal field 
distribution defined in Eq. ([3} for (a) g = 3, (b) g = 4. The numbers on each curve denote the value of A/J. The open 
circles represent the tricritical points, and the solid and dashed lines correspond to second- and first-order phase transitions, 
respectively. 
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FIG. 6; (Color online) (a) Temperature dependence of magnetization curves corresponding to Fig. ((S^) for g = 3 (a) with 
D I J — —1.0 and for some selected values of A/J, and (b) with A/J = 6.0 and for some selected values of D/J. 




FIG. 7: (Color online) (a) Temperature dependence of magnetization curves for q = 4 corresponding to Fig. ((5)d) for (a) 
D/J — 2.0 and (b) D/J — —4.0 with some selected values of A/J. 



19 



1 n 




\ 

\ 


- (a) 0.01 

- (b) 0.05 
-(c) 0.1 

- (d) 0.2 


0.8 


- 


0.6 


- 


1 




0.4 


■ .1 


n\ 




0.2 


c 
■ D/J=2.0 


/ 




n n 


.q=3 


L. 











1.0 



0.9 - 



0.8 - 



0.7 - 



0.6 - 



0.5 





D/J=2.0 


' ^ 


c 

V 


1 

V 


a q=4 

(a) 0.01 

(b)0.05 

(c)0.1 

i (d)0.2 


• 




-H=^ 






1 





3.0 3.5 4.0 4.5 5.0 5.5 6.0 



5.0 



5.5 



6.0 



6.5 



A/J 



A/J 



FIG. 8: (Color online) Variation of the magnetization curves as a function of A/J with D j J = 2.0. Left and right panels are 
plotted for q = 'i and 4, respectively. Each curve is plotted for different temperature values, namely ksT/J = 0.01,0.05,0.1 
and 0.2. 



